This article contains case-studies illustrating the benefits of implementing workflows with CoRC. Example code to get started can be found in the tutorial articles on model entitiy management, task management and model building.
library(tidyverse)
library(parallel)
library(CoRC)
# helper to run tasks in parallel on all cores
mapInParallel <- function(data, fun, ..., .prep = {}) {
cl <- makeCluster(detectCores())
clusterCall(cl = cl, fun = eval, .prep, env = .GlobalEnv)
result <- parLapplyLB(cl = cl, X = data, fun = as_mapper(fun), ..., chunk.size = 50)
stopCluster(cl)
result
}This example loads the Kummer2000 - Oscillations in Calcium Signalling model. The model has 3 species which oscillate. These oscialltions can be visualized as a trajectory through a 3D space. The example does this once in a deterministic and once in a stochatic fashion.
loadSBML("http://www.ebi.ac.uk/biomodels-main/download?mid=BIOMD0000000329")
#> # A COPASI model reference:
#> Model name: "Kummer2000 - Oscillations in Calcium Signalling"
#> Number of compartments: 1
#> Number of species: 3
#> Number of reactions: 8
# Run 24 sec (2 Periods)
setTimeCourseSettings(24, intervals = 10000)
timeseries <- list(
deterministic = runTimeCourse()$result,
stochastic = runTimeCourse(method = list(method = "directMethod", use_random_seed = T, random_seed = 1))$result
)
# Create the same plot for both timeseries
plots <- map(
timeseries,
plotly::plot_ly,
type = "scatter3d",
mode = "lines",
x = ~ `G-alpha`,
y = ~ activePLC,
z = ~ Calcium,
color = ~ Time
)
unloadModel()
plots$deterministicplots$stochasticThis implements an example from the Condor-COPASI paper. The example illustrates advantages of parallel processing.
# Run 1000 stochastic time series possibly in parallel
loadModel("https://raw.githubusercontent.com/copasi/condor-copasi/master/examples/stochastic_test.cps")
#> # A COPASI model reference:
#> Model name: "Kummer calcium model"
#> Number of compartments: 1
#> Number of species: 3
#> Number of reactions: 8
# timeseries <- 1:1000 %>% map(~ runTimeCourse()$result)
timeseries <-
# Defines parallel evaluation:
mapInParallel(
# perpare worker (.prep),
.prep = quote({
library(CoRC)
loadModel("https://raw.githubusercontent.com/copasi/condor-copasi/master/examples/stochastic_test.cps")
setTimeCourseSettings(method = list(method = "directMethod", use_random_seed = T))
}),
# iteration data (1000 random seeds) is given as .x to iteration formula,
1:1000,
# iteration code (formula ~)
~ runTimeCourse(method = list(random_seed = .x))$result
)
# Combine all results and reshape the data
plotdata <-
timeseries %>%
bind_rows() %>%
group_by(Time) %>%
# calculate mean and sd for all time points
summarise_all(funs(mean, sd)) %>%
# gather all values so the column "name" identifies "a_mean", "b_sd" etc.
gather("name", "value", -Time) %>%
# split up information on species (a,b,c) and type of value (mean, sd)
separate(name, c("species", "type"), "_") %>%
spread(type, value)
print(plotdata, n = 6)
#> # A tibble: 2,403 x 4
#> Time species mean sd
#> <dbl> <chr> <dbl> <dbl>
#> 1 0 a 8.00 0
#> 2 0 b 8.00 0
#> 3 0 c 8.00 0
#> 4 0.05 a 7.06 0.249
#> 5 0.05 b 8.12 0.117
#> 6 0.05 c 5.60 0.442
#> # ... with 2,397 more rows
plot <-
ggplot(data = plotdata, aes(x = Time, y = mean, group = species, tt_sd = sd)) +
geom_ribbon(aes(ymin = mean - sd, ymax = mean + sd, fill = species), alpha = 1 / 4) +
geom_line(aes(color = species)) +
guides(fill = "none") +
labs(
x = paste0("Time (", getTimeUnit(), ")"),
y = paste0("Concentration (", getQuantityUnit(), ")"),
color = "Species"
)
unloadModel()
plotly::ggplotly(plot, tooltip = c("group", "x", "y", "tt_sd"))This implements an example from the Mendes2009 paper on COPASI use cases.
loadSBML("https://www.ebi.ac.uk/biomodels-main/download?mid=BIOMD0000000068")
#> # A COPASI model reference:
#> Model name: "Curien2003_MetThr_synthesis"
#> Number of compartments: 1
#> Number of species: 9
#> Number of reactions: 3
setSteadyStateSettings(calculate_jacobian = FALSE, perform_stability_analysis = FALSE)
# Cartesian product of the input values
scan <- cross_df(
list(
cysteine = 0.3 * 10 ^ seq(0, 3, length.out = 6),
adomed = seq(0, 100, length.out = 51)
)
)
scan <-
scan %>%
mutate(
# Calculate steady state fluxes for every row
ss_fluxes = pmap(., function(cysteine, adomed) {
setSpecies("Cysteine", initial_concentration = cysteine)
setSpecies("adenosyl", initial_concentration = adomed)
ss <- runSteadyState()
stopifnot(ss$result == "found")
ss$reactions$flux
}),
# The second flux is CGS
CGS = map_dbl(ss_fluxes, 2),
# The third flux is TS
TS = map_dbl(ss_fluxes, 3)
)
plot <-
ggplot(data = scan, aes(x = adomed, group = cysteine)) +
geom_point(aes(y = CGS, color = "CGS")) +
geom_point(aes(y = TS, color = "TS")) +
geom_line(aes(y = CGS, color = "CGS")) +
geom_line(aes(y = TS, color = "TS")) +
labs(
x = paste0("Adomed (", getQuantityUnit(), ")"),
y = paste0("Flux (", getQuantityUnit(), " / ", getTimeUnit(), ")"),
color = "Species"
)
unloadModel()
plotly::ggplotly(plot)This implements an example from the Mendes2009 paper on COPASI use cases. It is in many ways similar to the previous example but is written to run parallelized.
loadSBML("https://www.ebi.ac.uk/biomodels-main/download?mid=BIOMD0000000068")
#> # A COPASI model reference:
#> Model name: "Curien2003_MetThr_synthesis"
#> Number of compartments: 1
#> Number of species: 9
#> Number of reactions: 3
# 10000 repeats of steady state task with random cysteine and adomed
scan <-
# Defines parallel evaluation:
mapInParallel(
# perpare worker (.prep),
.prep = quote({
library(CoRC)
loadSBML("https://www.ebi.ac.uk/biomodels-main/download?mid=BIOMD0000000068")
setSteadyStateSettings(calculate_jacobian = FALSE, perform_stability_analysis = FALSE)
}),
# iteration data (10000 random seeds) is given as .x to iteration formula,
1:10000,
# iteration code (formula ~)
~ {
set.seed(.x)
cysteine <- 0.3 * 10 ^ runif(1L, min = 0, max = 3)
adomed <- runif(1L, min = 0, max = 100)
setSpecies(
key = c("Cysteine", "adenosyl"),
initial_concentration = c(cysteine, adomed)
)
ss <- runSteadyState()
stopifnot(ss$result == "found")
list(
cysteine = cysteine,
adomed = adomed,
CGS = ss$reactions$flux[2],
TS = ss$reactions$flux[3]
)
}
)
# Combine all results and reshape the data
plotdata <-
scan %>%
bind_rows() %>%
gather("reaction", "flux", CGS, TS)
plot <-
ggplot(data = plotdata, aes(x = adomed, y = flux, group = reaction, tt_cys = cysteine)) +
geom_point(aes(color = reaction), alpha = 1 / 10, size = 3 / 4) +
labs(
x = paste0("Adomed (", getQuantityUnit(), ")"),
y = paste0("Flux (", getQuantityUnit(), " / ", getTimeUnit(), ")"),
color = "Species"
)
unloadModel()
plotly::ggplotly(plot, tooltip = c("tt_cys", "x", "y"))This implements an example from the Mendes2009 paper on COPASI use cases.
loadSBML("http://www.ebi.ac.uk/biomodels-main/download?mid=BIOMD0000000010")
#> # A COPASI model reference:
#> Model name: "Kholodenko2000 - Ultrasensitivity and negative feedback bring oscillations in MAPK cascade"
#> Number of compartments: 1
#> Number of species: 8
#> Number of reactions: 10
# get timeseries for the record
data_before <-
runTimeCourse(1000, 1)$result %>%
set_tidy_names(TRUE)
# read experimental data
data_experimental <-
read_tsv("data/MAPKdata.txt") %>%
rename(Time = time, "Mos-P" = "MAPKKK-P", "Erk2-P" = "MAPK-P") %>%
set_tidy_names(TRUE)
# define the experiments for COPASI
fit_experiments <- defineExperiments(
data = data_experimental,
type = c("time", "dependent", "dependent"),
mapping = c(NA, "{[Mos-P]}", "{[Erk2-P]}"),
weight_method = "mean_square"
)
# define the parameters for COPASI
fit_parameters <-
map(
parameter_strict(regex(c("V1$", "V2$", "V5$", "V6$", "V9$", "V10$"))),
~ {
val <- getParameters(.x)$value
defineParameterEstimationParameter(parameter(.x, "Value"), start_value = val, lower_bound = val * 0.1, upper_bound = val * 1.9)
}
)
result <-
runParameterEstimation(
parameters = fit_parameters,
experiments = fit_experiments,
method = "LevenbergMarquardt",
update_model = TRUE
)
# get timeseries for the record
data_after <-
runTimeCourse(1000, 1)$result %>%
set_tidy_names(TRUE)
plots <- list(
Erk2.P =
ggplot(mapping = aes(x = Time, y = Erk2.P)) +
geom_point(data = data_experimental, aes(color = "experimental")) +
geom_line(data = data_before, aes(color = "before")) +
geom_line(data = data_after, aes(color = "after")) +
scale_color_manual(values = c(before = "red", after = "green", experimental = "black")) +
labs(
x = paste0("Time (", getTimeUnit(), ")"),
y = paste0("Erk2-P (", getQuantityUnit(), ")"),
color = "Series"
),
Mos.P =
ggplot(mapping = aes(x = Time, y = Mos.P)) +
geom_point(data = data_experimental, aes(color = "experimental")) +
geom_line(data = data_before, aes(color = "before")) +
geom_line(data = data_after, aes(color = "after")) +
scale_color_manual(values = c(before = "red", after = "green", experimental = "black")) +
labs(
x = paste0("Time (", getTimeUnit(), ")"),
y = paste0("Mos-P (", getQuantityUnit(), ")"),
color = "Series"
)
)
unloadModel()
result$fitted_values
#> # A tibble: 2 x 5
#> fitted_value objective_value root_mean_square error_mean error_mean_std…
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 [Mos-P] 25.6 1.69 -0.287 1.66
#> 2 [Erk2-P] 10.0 1.05 0.339 0.999
result$parameters
#> # A tibble: 6 x 8
#> parameter lower_bound start_value value upper_bound std_deviation
#> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 (MAPKKK … 0.25 2.47 2.47 4.75 0.140
#> 2 (MAPKKK … 0.025 0.247 0.247 0.475 0.00831
#> 3 (dephosp… 0.075 0.886 0.886 1.42 0.260
#> 4 (dephosp… 0.075 1.42 1.42 1.42 0.708
#> 5 (dephosp… 0.05 0.720 0.720 0.95 0.0799
#> 6 (dephosp… 0.05 0.698 0.698 0.95 0.0956
#> # ... with 2 more variables: coeff_of_variation <dbl>, gradient <dbl>
result$protocol
#> [1] "Algorithm started at 2018-07-31 22:13:12.\nFor more information about this method see: http://www.copasi.org/tiki-index.php?page=OD.Levenberg.Marquardt\n\nIteration 164: Objective function value and parameter change lower than tolerance (1/3). Resetting lambda.\n\nIteration 165: Objective function value and parameter change lower than tolerance (2/3). Resetting lambda.\n\nIteration 166: Objective function value and parameter change lower than tolerance (3/3). Terminating.\n\nAlgorithm reached the edge of the parameter domain 299 times.\n\nAlgorithm finished at 2018-07-31 22:13:13.\nTerminated after 167 of 2000 iterations.\n\n"
plotly::ggplotly(plots$Erk2.P)plotly::ggplotly(plots$Mos.P)